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(54) Efficient cone-beam reconstruction system using circle-and-line orbit data 

(57) A system (10) is provided for reconstructing x- 
ray tomographic images based upon the circle-and-line 
orbit. The system includes the steps of combining, in the 
frequency domain, images reconstructed separately 
from the circular and linear orbit data. Specifically, the 
method involves the steps of (1) converting the shadow 
zone (28) into a shadow cone (30) in frequency domain, 
(2) mapping the shadow cone in 3D Fourier space onto 
a set of 2D Fourier planes (27), (3) removing data lying 
within the region marked as the projection of the 
shadow cone from the Fourier transform of the circular 
orbit data and patching the same from the Fourier trans- 
form of the line orbit data, (4) transforming each 2D Fou- 
rier plane into a respective 2D image plane, and (5) 
converting the horn-shaped volume back to a grid vol- 
ume. An interpolation technique is also provided for 
reconstructing the line orbit data using a Direct Fourier 
method (DFM). 
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Description 

[0001] The field of the invention relates to computer tomography and more particularly to x-ray medical computer 
tomography. 

5 [0002] Computer tomography (CT) using cone-beam projection technology typically employs an x-ray source to 
form a vertex in the shape of a cone. A scanning orbit of the x-ray source is a curve along which the vertex of the cone 
beam moves during a scan. A set of detectors are disposed at a fixed distance from the cone-beam source. Methods 
currently exist which have established a relationship between cone-beam projections and a first derivative of the three 
dimensional (3D) Radon transforms of such projections. Theoretically, each source position "S" can deliver Radon data 

w positioned on a sphere, the Radon shell of "S." The corresponding spheres on a circular orbit sweep out the Radon 
space of the object, referred to as a torus. Regions inside the smallest sphere containing the object, but exterior to the 
torus, is called the shadow zone. The shadow zone characterize the region where Radon data is missing from the cir- 
cular orbit. 

[0003] When using a circle-and-line orbit, the line-orbit data is typically used for filling in the shadow zone. Even 
15 though the shadow zone is very small, as compared to the volume of the entire valid Radon data delivered by the cir- 
cular orbit, current cone-beam projection techniques are computationally intensive. Accordingly, a need exists for sim- 
plifying image reconstruction based on the circle-and-line orbit. 

[0004] According to the present invention, a system is provided for reconstructing x-ray tomographic images based 
upon the circle-and-line orbit. The system includes the steps of combining, in the frequency domain, images recon- 

20 structed separately from the circular and linear orbit data. Specifically, the method involves the steps of (1) converting 
the shadow zone into a shadow cone in the frequency domain, (2) mapping the shadow cone in 3D Fourier space onto 
a set of 2D Fourier planes, (3) removing data lying within the region-marked as the projection of the shadow cone from 
the Fourier transform of the circular orbit data and patching the same from the Fourier transform of the line orbit data, 
(4) transforming each 2D Fourier plane into a respective 2D image plane, and (5) converting the horn-shaped volume 

25 back to a grid volume. An interpolation technique is also provided for reconstructing the line orbit data using a Direct 
Fourier Method (DFM). 

[0005] An embodiment of the invention will now be described, by way of example, with reference to the accompa- 
nying drawings, in which: 

30 FIG. 1 depicts an x-ray tomography imaging system of the present invention; 

FIG. 2 depicts a cross-section of the circular orbit formed by the system of FIG. 1 ; 

FIG. 3 depicts a shadow zone of FIG. 2 and corresponding shadow cone; 

FIGs. 4a and 4b illustrate a cross-section of the shadow cone of FIG. 2 projected onto a plane; 

FIG. 5 depicts incremental processing steps performed in the z-direction by the system of FIG. 1 ; 
35 FIG. 6 depicts a set of parallel projection data at the 0 max angle (obtained from the M th row of the detector array as 

depicted in FIG. 5); 

FIG. 7 depicts a set of parallel projection data at the 0 min angle (obtained from the first row of the detector array as 
depicted in FIG. 5); and 

FIG. 8 depicts inverse bilinear interpolation of the system of FIG. 1 . 

40 

[0006] FIG. 1 depicts a CT imaging system 10, generally, in accordance with an illustrated embodiment of the 
invention. Included within the system 10 is a central processing unit (CPU) 16, a sampling device 12, and other image 
support and processing subsystems 14, 18, 20, 22. The sampling device 12 includes an x-ray source 34 and an image 
detecting array 32 along with actuating devices. The image support devices 14, 18, 22 include a Fast Fourier processor 

45 14, an interpolator 18, a memory 20, and a display 22. 

[0007] In two dimensional (2D) tomographic reconstruction one alternative to Filtered Back Projection (FBP) is 
Direct Fourier Method (DFM) back projection. Direct Fourier Method back projection is based on the projection slice the- 
orem, which postulates that a one dimensional (1D) Fourier Transform (FT) of a projection at a specific angle corre- 
sponds to the cross section of the 2D FT of the object at the same angle. DFM shall be used as an alternative to the 

so conventional back-projection reconstruction involving the line orbit. 

[0008] In a CT system 10, (FIG. 1) using circular orbit data and alternatively the circle and line orbit data which is 
to be reconstructed, represented in the function defined by r"(7), can be expressed by the following equation: 

55 

[0009] The first two elements, f cQ (r) and f c1 (7), are evaluated using circularly-scanned data and techniques (e.g., as 
given by Feldkamp et al. and Hu). For example, Feldkamp et al. has published cone-beam algorithms in a paper enti- 
tled, PRACTICAL CONE-BEAM ALGORITHM, published in the Journal of the Optical Society of America, Vol. 1 , No. 6, 
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June 1984, pages 612 - 619, herein incorporated by reference. The second element is needed only when the third term 
is evaluated In the Radon space. Since we will evaluate the third term in the Fourier space, thus the second term will 
be ignored. 

[0010] The third element, f}(7), is evaluated from the linear orbit data. It is typically used for filling in a shadow zone 
5 28 defined by a torus 24, which is swept out by the circular orbit about the axis of rotation, and the smallest sphere con- 
taining the object 40 (FIG. 2). Even though this shadow zone 28, which is defined in Radon space, is small as compared 
to the torus 24, the current back-projection process is computationally intensive, as is apparent from the complexity of 
f,(7) as defined, for example in a paper by Hui Hu, entitled, CONE BEAM RECONSTRUCTION ALGORITHM FOR THE 
CIRCULAR ORBIT AND A CIRCLE-AND-LINE ORBIT, published by the Applied Science Laboratory, GE Medical Sys- 
10 terns, 1/31/97, herein incorporated by reference, in equation 1 listed below. 

ji 

f ' Cr) = 2 J ~ ~, ldz o J deH ( z o.©.0,^ oSi ne + z 0 cose. 0> 

47t 2 (d+r-x l ) e J =Q 



15 



25 



30 



where 



20 H(z o ,0,/) = cos0*v(z o ,0,/) 



^ 2 + / 252z z 0 (/.e) +2/ 6i Zo (/,e)>| 



d 2 5 2 / d 2 5 ' ) 



£ ZQ (l,e) = JJ P Zo ( Y,Z)8( Vsine+ZcosG-/) dY dZ , 



/ - .\ f I when2Iz a cos0 + z\ cos 2 &-d 2 sin 2 0 > 0 
[0 otherwise 



and where H(z o ,0,/) is the weighted second derivative of Radon data computed from a snapshot incurred at 
z = z 0 . For each point in the object space, r f the point is associated with a point ( Y^Zq) on the detector plane. If a circle 
35 is drawn on the detector plane using the origin (0, 0) and (V^,Z 0 ) as two ends of a diameter, then all the valid Radon 
data, H(z 0 , 0, I), intercepted by the circle satisfy the constraint: /=V o sin0+Z o cos0 . For each point (Y&Zq), there is a 
corresponding set of Radon data, /-/(z o ,0,/), to be summed and then back-projected onto the object space. The com- 
putational complexity of this process is thus close to 0(N 5 ). 

[0011] The present invention describes a method for determining a computational efficient CT reconstruction using 
40 data obtained from a line orbit scan of the object 40 (illustrated in FIG. 2). First, the shadow zone 28 in Radon space is 
augmented and converted into a shadow cone 30 in frequency space (FIG. 3). Secondly, the patching process per- 
formed by CPU 16 is simplified from the original three-dimensional frequency space into multiple 2D frequency slices. 
[0012] The first step is based on the fact that the Fourier transform of a radial line in 3D Radon space is equivalent 
to the same radial line in 3D Fourier space. If a radial line in 3D Radon space intercepts the shadow zone 28, then its 
45 counterpart in the 3D Fourier space is "contaminated" because some Radon data in the radial Radon line has been 
missing. Since all the radial lines intercepting the shadow zone 28 form a shadow cone 30, the collection of their coun- 
terparts in the 3D Fourier space also form a frequency-insufficient cone. The frequency-insufficient cone (or the con- 
taminated cone) of the first term of equation (1), which is defined in Fourier space, will be replaced by that of the third 
term of equation (1). 

so [001 3] The second step is to show that the linear orbit supports enough frequency data to replace the contaminated 
cone (which is due to the circular orbit), followed by a procedure for replacing the contaminated cone in 2D Fourier 
space. 

[0014] FIG. 5 shows that the linear scanning orbit may be implemented by holding an x-ray source 34 and a detec- 
tor panel 32 stationary, while moving object 40 along the z-axis 36 (which is also the rotational axis of the circular orbit). 
55 The object 40 is shown to move in the M z direction" in discrete increments identified at u z_step" 36 increments. Each 
vertical column of detectors 32 together with the source 34 defines a plane through the object 40. Within each vertical 
plane, every detector position defines a specific angle and the data collected by that detector corresponds to a parallel 
projection of the object sampled by intervals of z-step 36. FIG. 6 shows a number of parallel projections of object 40 at 
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angle 9 m =9 max and FIG. 7 shows parallel projections of object 40 at angle e m =0 min . These parallel projections are 
used to fill in the corresponding missing angular segment in the 2D FT plane. 0 max and 0 min are angles selected so that 
the data volume for the shadow cone 30 (FIG. 3) used to replace the shadow zone 28 data volume encompasses the 
shadow zone data volume. This condition is satisfied when G o =0 m , as illustrated in FIG 2, where 

e 0 = sin- 1 (£), 



10 wherein "R* is the radius of the smallest sphere containing object 40 and "d" is the radius of the circular orbit. 

[0015] The shadow cone 30 (FIG. 3) in 3D Fourier space may be projected as a 2D shadow cone onto a set of 2D 
Fourier planes. The inverse transform of each of these 2D Fourier planes corresponds to an image plane defined by the 
x-ray source and a column n N" of detectors 32 (FIG. 5). 

[0016] First, let F(u,v) be the Fourier transform of a 2D function f(x,y). The 1 D Fourier data on x = x 0 is derived as 
15 follows: 

F Xq (v) = jF(u,v)e i2nXQV du = Ja/^(x,y)e- y2n(xu ^ ) dxdy}e /27CJf0 "c/u (2) 

= /Jf(x,y)e^ 2 ^{Je^ 27c(x - Xo)u du}dxdy 

= M*,y)8(x-x 0 )e ' J2ny * dxdy = Jf(x 0 , y)e J2nyv dy. 

[0017] Similarly, for the Fourier transform of a 3D function f(x, y, z), the 2D Fourier transform of the object slice at 
25 x = x Q is identified as follows: 

F Xq (v,w) = jF(u,v, W )e i2nX ° V du = J{//f(x,y,z)e- y2 ^ (xu+ ^ +zw) dxdyc/z}e /27IXot; du (3) 
= J/f(x,y,z)e- y2lc(yv+zlv) {Ie" y27l(;f " Xo)U du}dxdydz 
=//f(x,y,z)8(x-x 0 )e- y2 ^ +zw) dxdydz = l^y^e^^dydz. 

[0018] For the case when the desired "n" dimensional Fourier space is not orthogonal to the (n+1) th dimension of 
the augmented (n+1) th dimensional Fourier space. For example, if the desired 1D Fourier transform lies on a line 
35 defined by the equation xcos0+ysin0-p=O, it is appropriate to rotate the original 2D coordinates, (x,y), into new coor- 
dinates, (x',y'), by angle 0. Let F(u' v') be the Fourier transform of f(x',y'). Then F(u',v') is also obtained by rotating F(u,v) 
by the angle 0. The Fourier transform along the line, xcos0+ysin0-p=O , is identified as follows: 

F x , p (v') = iF(u',v')e J2npul du\ (4) 

40 

[0019] Similarly, the 2D Fourier transform along the plane defined by x cos cp + y sin <p - p = 0 is given by 

F x .. p (v',w) = jF(u'y,w)e J2npu 'du\ (5) 

45 [0020] Again, F(u'v'w) is obtained by rotating the (u,v) coordinates by an angle q>, as illustrated in FIG. 4a, while 
keeping the w-coordinate unchanged. The rotated coordinate system is graphically illustrated in FIGs. 4a and 4b, where 
the new coordinates, (u',v\w), are chosen such that u' is parallel to the normal of the plane of interest (shown as the 
SO* line). As such, the shadow disc 26 (which is a cross-section of shadow cone 30 (FIG. 4a), intersecting a horizontal 
plane 27 (FIG. 4b), is projected onto SO' as a line segment centered at "O"' and having end points T and u e. u The 

so length of the line segment le equals the diameter of the shadow cone 30 centered at "O." 

[0021] In summary, shadow cone 30 (FIG. 4a), defining the area in Radon space in 3D Fourier space, is projected 
onto 2D Fourier subspace by applying a 1D inverse Fourier transform. In other words, rather than by replacing the 
shadow cone in 3D Fourier space, the scanned object 40 is represented in a series of 2D Fourier planes, and replaces 
the projected 2D shadow disc 26 in all the individual planes. 

55 [0022] In applying the principle described above, the following array and scanner parameters will be assumed. 
Detector array 32 (Figure 5) will be referred to by the expression, D(M x N), where the detector array 32 has "M° rows 
and H N° columns. Projection data of circular orbit will be referred to by the expression P C (M x N x V c ), where V c is the 
total number of views. Projection data of the line orbit will be referred to by the expression P,(M x N x V^), where V, is 
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10 



also the total number of views. An image reconstructed from circular orbit data will have the form, l c (K x K x K) or l C s( K 
x N x K); the former is a cubic volume, while the latter is a "horn-shaped" volume defined by the x-ray source and the 
detector panel. An image reconstructed from line orbit data will have the form, l|(K x K x K) or l| S (K x N x K); likewise, 
the former is a cubic volume and the latter is a "horn-shaped" volume. 

[0023] The image is reconstructed from the line-orbit data. As a first step, the line-orbit projection data is modified 
to accommodate the effect of divergence as a plurality of rays passes from source 34 to detector array 32, as graphi- 
cally illustrated in FIGs. 5, 6, and 7. The divergent ray detected on the detector panel is modified by multiplying the 
detector values with a weighting function as follows: 



W{m,n) = 



jD 2 +(m-m 0 ) 2 a 2 +(n-n Q ) 2 b 2 



where n D a in this case is the source-to-detector distance (FIG. 5), and wherein ISO is projected onto (m 0 ,n 0 ) of 
15 the detector panel 32. The values "a" and "b" are the detector pitches in the row and column dimensions, respectively. 
Detector pitches are standard definitions defining the distance between the centers of two adjacent detectors. 
[0024] The n th angular vertical plane is reconstructed next. To reconstruct the n th angular vertical plane, the projec- 
tion data, P /n (M x Vj), is extracted from the n th column of the P/ array, where the m th row corresponds to a fixed beam 
angle of 0 m . 

20 [0025] The row elements of P m are shifted such that the object center of this vertical slice is realigned at 

25 

Let P in be the realigned projection array. 

[0026] To return to object space, either one of two methods may be taken. As a first option, apply a 1 D FFT to each 
row of Pin . Place each set of FFT data on a radial line at an angle 9 m for m = 1 ... M. An inverse FFT is utilized to con- 
vert the Fourier data set to object space. 
30 [0027] As a second option the filtered P\ n for all "N" columns are back-projected to obtain a horn-shape volume 
comprising "N" vertical planes defined by the x-ray source 34 and D N n detector columns "N." Let f n (x,y) be the recon- 
structed image based on the n th column, and Q n the filtered projection array, then 



M 



35 f n(*>y)=^M^I, O nm (xcosQ m + ys\nQ m ). (6) 

m=1 



[0028] The results of applying the above equation to the shifted row elements are identified as a volumetric array 
40 1 1s , having a horn shape. 

[0029] The final image may be reconstructed by one of two possible methods. First, the shadow cone 30 may be 
replaced in 3D Fourier space. Alternatively, the shadow cone 30 may be replaced in 2D Fourier spaces and then inter- 
polated from the horn shape volume back to cubic volume. 

[0030] The replacement of the shadow cone 30 in 3D space is presented first. As a first step, the Feldkamp algo- 
45 rithm is used to reconstruct the image, l c , associated with torus 24, from the circular orbit data. Next, a 3D Fourier trans- 
form, F(l c ), is obtained from l c . 

[0031] The horn shaped volume, l| S , is interpolated to obtain a corresponding cubic volume, l ( . A 3D Fourier trans- 
form, F(l|) is obtained from I,. 

[0032] The shadow cone portion of F(l c ) is replaced with its counterpart from F(l|) to provide a repaired Fourier data 
so array. An inverse Fourier transform is then applied to the repaired Fourier data array to return to cubic space. 

[0033] Replacing the shadow cone in 2D Fourier spaces is presented next. As a first step, the Feldkamp algorithm 
is again used to reconstruct image l c from the circular orbit data. Next, l c is interpolated to l cs , in which l s defines a plu- 
rality of angular vertical planes. I| S is obtained from the line orbit data, as described above. 

[0034] The Fourier transform F(l cs ) and F(l iS ) is obtained from l cs and l )S - The shadow cone is substituted into each 
55 2D plane of F(l cs ), with its counterpart defined in F(I )S ). An inverse Fourier transform is applied to each of the repaired 
2D sets of Fourier data. Finally, the horn-shape volume is interpolated back to cubic volume. 

[0035] The DFM method described above requires interpolation from polar coordinates to rectangular coordinates. 
[0036] The cone beam reconstruction described above is based on circle and line orbits which requires that the 
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45 



missing shadow zone 28 of circle orbit data in the Fourier domain to be replaced by the line orbit data. For the line orbit 
orientation illustrated in FIG. 5, the object 40 is shown to move in the "z direction" in discrete increments identified as 
"z_step" 36 increments. Each vertical column of detectors 32 together with the source 34 defines a plane through the 
object 40. Within each vertical plane, every detector position defines a specific angle and the data collected by that 

5 detector corresponds to a parallel projection of the object sampled by intervals of z-step (FIGs. 5, 6, and 7). FIG. 6 
shows a number of parallel projections of object 40 at angle 0 m =0 max and FIG. 7 shows parallel projections of object 
40 at angle e m =e min . These parallel projections are used to fill in the corresponding missing angular segment in the 
2D FT plane. 0 max and 0 m | n are angles selected so that the data volume for the shadow cone 30 (FIG. 3) used to replace 
the shadow zone 28 data volume encompasses the shadow zone data volume. This condition is satisfied when 

10 0o =0 m , as illustrated in FIG 2. 

[0037] One approach is to reconstruct the filtered back-projection (FBP) image for each angular plane using the 
parallel projections from the set of limited angular views and then to compute the 2D FT and use it to replace the miss- 
ing angular segment. 

[0038] Another alternative is to convert the polar grid data from 1D FT of parallel projections to rectangular grid 
15 directly in 2D FT domain. This method is similar to the interpolation required in Direct Fourier Method back-projection. 
[0039] In the line orbit, for each plane, a number of projections (equal to °M" detectors in each vertical column) are 
acquired, for an angular range from 0 mjn to 0 max> Interpolation techniques for the line data provide the methodology for 
developing the below described interpolation techniques that are useful for calculating the cone beam reconstruction 
methods. 

20 [0040] First, let ^(x, y) be the image in a fixed coordinate system x - y and u^x , y ) represent the object in a coor- 
dinate system x - y rotated from x - y by an angle <t>. The projection of the object at a view angle <t> is 

P^(x) = i[i^x f y)dy. (7) 

25 [0041] The 1 D FT of p^(x ) is given as 

p*(u) = Jp*(*y 32n * u du. (8) 

30 

[0042] If Ivyp, <|>) denotes the 2D FT of U-^x, y) in polar coordinates, then based on the Projection slice theorem, 

M(p, <t>) = P^(p), 

35 

where 

p > 0, 0 < <t>< 2k (9) 

40 [0043] Assuming that M(p, (j>) is known, this FT data on polar coordinates may be converted to Cartesian coordi- 
nates by the expression, M(p cos p sin <|>) = M(p, <t>) . The object is reconstructed by 2D inverse FT using the 
expression, 



H(x, y) = f M c (u, v)e - j2{ ^ x+vy) dudy. (1 0) 



[0044] However, M(p, is not known at all positions but only at finite set of discrete points (pj, <|) k ). The problem 
then becomes interpolating from the known values of polar points to the values requires over a rectangular Cartesian 
grid. 

[0045] An existing theorem is utilized to solve this problem. First, let x(t) be a periodic function in time space V with 
so period "T." If x(t) is angularly band-limited to a value (c| ) where "c M is a constant, and "K" is an integer relating to the 
number of samples of x(t), then x(t) can be written as: 



55 
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1-0 



10 T 



*[f (t-lli)] 
x(t) = ^(iT,) L J 



(2K + l)sin[|(t - IX,)]' 



where T s = . (11) 

s 2K + 1 



15 

[0046] The proof follows from Shannon's sampling theorem as follows: 

*(t) = t xta>inc£i (t - n%)j (12) 



20 



25 



30 



35 



40 



= - + g x((l - N)T K >in^ (t - (1 - N>Tjj 
+ g x(lX,>incfi (t-lTjl 

1-0 L x s J 

+ X *((1 + N)T s >in^ (t - (l + N)T x )j + ... 

= j£ x(lT s >inc^i (t - (IX, + kT))j (13) 



45 were the number of data array columns N = 2K + 1 . 

[0047] Angular interpolation employs the function M(p, wherein variable $ has a period of 2n. Since M(p, <j>) is 
angularly band-limited to n K B (i.e., Fourier series coefficients satisfy c n = 0 for Inl > K), then 



2K+1 
k=0 



where 

55 



m(p.«- X m(p^X*"I&) (14) 
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35 



(2K + 1) si 



(15) 



10 



[0048] Radial interpolation is described next. If the object is space limited to diameter 2A, where "A B is equivalent 
to radius "R", identified previously as the radius of the smallest sphere containing the object 40, its projections are also 
limited to diameter 2A and therefore the FT of the projection P<|>(u) can be reconstructed as 



15 



■v«>- i ^HK"- Ai- 



de) 



20 using Shannon's sampling theorem. 

[0049] Combining angular and radial interpolation: 



25 



2K+1 



(17) 



n=-<» k=0 



[0050] M C (U|, v m ) is computed over the Cartesian space as: 



30 



M c (u,, vj = M^ul + v;, cos" 1 U " 



(18) 



40 



[0051] To reduce interpolation error, a number of steps may be taken, as described below. Error from inadequate 
sampling in space will be considered first. 

[0052] During the line orbit, the parallel projection data are sampled in intervals of a predetermined z-step 36. As 
the projections are space limited they will not be band-limited. Sampling the projections over a predetermined z-step 
window (for example 2A) will effectively reduce the bandwidth and introduce aliasing artifacts into the 2D FT domain. 
Consequently, the n z_step u 36 is selected to be small enough such that the effective bandwidth, 



45 



1 



2*z_step 



[points/cycle], 



will substantially include the frequency content of the projections. 

[0053] Another source of error is inadequate sampling in a radial frequency. If the object is space limited in diameter 
so to ^ , its projections will also be space limited to 2A, as such, the FT is sampled uniformly at intervals of at least, ^ 
[points/cycle], in order to avoid aliasing. 

[0054] For the line data, the effective bandwidth is defined as "B." The number of radial samples "N" are chosen 
such that 



55 



2TB _ __L_<J_ 
N N*z_step~2*A 
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Also the grid spacing in the 2D FT domain is at most ^ [points/cycle]. 

[0055] Another source of error is an insufficient number of views. Non-isotropic and off-centered objects will not be 
angularly band-limited. So, the number of projections may be very large. However, in practice the number of projections 
is chosen such that in the 2D FT domain the sampling between two neighboring radial lines will be approximately equal 

5 to the maximum grid spacing of ^ [points/cycle]. 

[0056] Interpolation in the 2D FT domain is described next. For the line orbit data, a rectangular grid spacing of ^ 
[points/cycle], is used in the 2D FT domain. In a computer implementation, typically, the Fast Fourier Transform (FFT) 
is used for computing the 1 D FT of each projection. In order to improve radial sampling, zero padding of the projection 
data is performed before applying FFT. In this specification zero padding involves the steps of adding zeros to the data 

w set to minimize artifacts after an FFT is performed. Once a radial sampling has been obtained, Nearest neighbor inter- 
polation or Linear interpolation is used to compute the values between the new sample points. 

[0057] As mentioned previously, for the line orbit a plurality of projections are required for a small angular segment 
(i.e., angular direction is over-sampled compared to radial sampling). One concern is how to use this data efficiently A 
simple interpolation scheme such as nearest neighbor will not use all of the available data. One alternative approach is 

15 to decimate the number of views by grouping the rows and then summing the rows from each group so as to increase 
the signal to noise ratio (SNR) utilizing an inverse bilinear interpolation technique, as further described below. 
[0058] In the present invention, an inverse bilinear interpolation technique is used on data in the polar grid format. 
For optimal interpolation in the azimuthal direction, a(<(>) the function defined in equation (13) may be used. However, 
this function is computationally impractical. For a computer implementation, values in the interpolation function must be 

20 truncated which causes artifacts. 

[0059] From a computational standpoint, the simplest interpolation algorithm to implement is the Nearest neighbor 
(NN) algorithm, wherein each pixel is given the value of the sample which is closest to the pixel. NN is equivalent to con- 
volving the sampled data by a rectangular function. Convolution with a rectangle function in the spatial domain is equiv- 
alent to multiplying the signal in frequency domain by a sinc(. ) function, which is a poor low pass filter since it has 

25 prominent side lobes. 

[0060] Linear interpolation amounts to convolution of the sampled data by a triangle window. This function corre- 
sponds in frequency domain to a sinc 2 (. ) function which has a better low pass filtering property than NN. 
[0061] The NN method interpolates on the basis of a single point. The linear interpolation algorithm interpolates on 
the basis of two nearest points. Using three points for interpolation would result in two points on one side and only one 
30 on the other side. The B-spline method enables interpolation over the two nearest points in each direction. The Cubic 
B-spline method comprises four convolutions of the rectangular function of the NN method. Since Cubic B-splines are 
symmetric, they only need to be defined in the interval (0, 2). Mathematically, the B-splines can be written as: 

f(x) = x 3 -x 2 + | xe(0,1) 



f(x) = ^+x 2 -2x+| xe (1,2) 

40 

[0062] The Cubic B-splines equations defined above may be adapted to fit a smooth third-order polynomial through 
the available set of shadow cone 30 data. 

[0063] Another interpolation function involves the truncation of a(<|>). As mentioned previously, it is desirable to find 
the ideal interpolation function in the angular direction of the cs(ty) function. However, the a(<|>) function has considerable 
45 energy and, as such, does not decay fast over an extended distance, and therefore cannot be easily implemented as a 
space domain convolution. It is natural to truncate the function over a small distance, but truncation disregards much of 
the above said energy. Moreover, truncation in the space domain will produce ringing artifacts when the object is con- 
verted from the space domain to the signal domain. 

[0064] Bilinear interpolation is a 2D interpolation method. In this interpolation method f(x,y) is evaluated by a linear 
so combination of f(n 1f n2) at the four closest points for n 1 T < x < (n 1 + 1)T and n 2 T < y < (n 2 + 1)T . The interpolated 
value f(x,y) in the bilinear method is: 

f(x, y) = (1 -A x )(1 - A y )f(n 1f n 2 ) + (1 - A x )A y f(n 1 , n 2 + 1) + A x (1 -A y )f(n 1 + 1, n 2 ) + A x A y f(n 1 +1,n 2 + 1) 

55 where 
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(x- ni T) 



and 



(y-n 2 T) 



10 



15 



[0065] Inverse bilinear interpolation converts polar-grid data into rectangular-grid data via inverse bilinear weight- 
ing. That is, each polar-grid element is redistributed to its four nearest rectangular grids using the four sub-areas as nor- 
malized weights. The method is robust in the sense that the original data does not have to be regularly distributed, as 
long as the data has at least one element in each rectangular square. As such, undesirable data elements caused by 
defects or seams in detector array 32 (FIG. 5) may be ignored. Another benefit of this process is that ail the good data 
are fully utilized, which, consequently, provides a better signal-to-noise ratio (SNR). 

[0066] A linear orbit datum value H x B and associated grid points used in inverse bilinear interpolation are graphically 
illustrated in FIG. 8. The data value "x" is distributed to its four nearest grid points as follows: 



20 



25 



30 



35 



40 



45 



50 



9ij = 9ij + au> 



9jj+\ - 9jj+\ + a 10 x 



9i+\j- 9j+\j+ a 01 x 
9j+\j+\ = 9j+\j+\ + a oo x 



[0067] Each grid point has an accumulated value, gy, which identifies the accumulated weighted °x 8 values from the 
four adjacent squares. To normalize the accumulated value at each grid point, it is necessary to maintain another vari- 
able, Wjj, which identified the accumulated weights (a,y) associated with the bilinear distribution of pixel (/J). In fact, this 
is a normalization step at the end of the interpolation procedure because each grid value is divided by accumulated 
weight, Wjj. 

[0068] Six sets of simulations were conducted using the described interpolation methods for the 2D FT domain and 
the image domain. The simulation results illustrate Inverse bilinear interpolation as compared to Nearest neighbor (NN), 
Linear, Cubic B-Spline, and Filtered back projection (FBP) interpolation methods. Results are presented in Tables I 
through VI. 

[0069] For the first three simulations, Linear interpolation was used in the radial direction. Nearest neighbor, Linear 
and Cubic B-spline interpolation techniques were used in the azimuthal direction. The reconstruction of an angular seg- 
ment in the 2D FT domain was simulated using the line orbit data. The slice of interest was the central slice through the 
object in FIG. 5. The object consisted of three horizontal disks so that it had sharp edges in z (vertical) direction. 
[0070] The Truncated o(§) interpolation method was also used in the azimuthal direction (including 8 neighboring 
points). Linear interpolation was used in the radial direction. Nearest neighbor interpolation and Linear interpolation 
generated better performance results than Cubic B-spline interpolation. However, Cubic B-spline uses only 4 points as 
opposed to 8 points in Truncated o((J>) interpolation. So, Cubic B-spline is computationally more attractive than the 8 
points Truncated a(<|>) interpolation, as such, Truncated a((j>) interpolation was not included in the comparison. 
[0071 ] Polar points based bilinear interpolation was also evaluated. The 2D FT of the slice of interest was computed 
and used as the base reference for of comparing the above identified methods. The 2D FT of the FBP reconstruction 
was also included in the comparison. FBP reconstruction requires the reconstruction of the 2D FMP image and then 
computation of the 2D FT. However, Polar bilinear interpolation does not require conversion back to the image domain. 
This results in a computational advantage over the FBP method. 

[0072] In order to measure the performance of each method, the interpolated rectangular grid values were com- 
puted using the above methods for different numbers of views and radial samples. The error was computed with the true 
2D FT as follows: Norml (total absolute difference, Table I); Norm2 (square root of the total squared distance, Table II); 
and Normoo (maximum absolute difference, Table III). Note, that terms Norm 1 , Norm2, and Norm°o are standard terms 
used in statistical analysis. 
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Table I 





% NORM1 ERROR in 2D FT DOMAIN 


# of views 


256 


256 


64 


16 
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Table I (continued) 
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% NORM1 ERROR in 2D FT DOMAIN 


n ot raaiai samples 


1 tou 








Mpor Mpinh 

INCCtl. INCiyil. 


58.4 


58.0 


60.9 


75.2 


Linear 


49.7 


49.5 


51.6 


67.1 


Cubic B-Spline 


45.7 


45.6 


47.9 


64.9 


Inverse bilinear 


41.2 


40.7 


43.3 


60.9 


FBP 


77.6 


77.6 


77.6 


77.6 



15 

Table II 



25 





% NORM2 ERROR in 2D FT DOMAIN 


# of views 


256 


256 


64 


16 


# of radial samples 


1280 


256 


256 


256 


Near. Neigh. 


31.6 


31.6 


33.0 


38.9 


Linear 


26.5 


26.5 


26.3 


33.8 


Cubic B-Spline 


24.4 


24.5 


24.8 


32.9 


Inverse bilinear 


21.8 


21.3 


22.4 


31 .9 


FBP 


54.5 


54.5 


54.5 


54.5 



Table III 





% NORMoo ERROR in 2D FT DOMAIN 


# of views 


256 


256 


64 


16 


# of radial samples 


1280 


256 


256 


256 


Near. Neigh. 


8.88 


8.88 


8.88 


8.88 


Linear 


8.87 


8.87 


7.82 


6.78 


Cubic B-Spline 


8.61 


8.61 


7.70 


6.77 


Inverse bilinear 


7.58 


6.08 


6.08 


7.26 


FBP 


80.25 


80.25 


80.25 


80.25 



[0073] The corresponding limited view images were also computed (Tables IV-VI) using inverse 2D FT from the 
interpolated rectangular grids. The errors were compared in the reconstructed images with the true image. The FBP 
so image are also included in the comparison. Note that in the real implementation for the line orbit, it is not necessary to 
reconstruct the images before patching the 2D FT data into the Feldkamp 2D FT data. 



Table IV 
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% NORM1 ERROR in IMAGE DOMAIN 


# of views 


256 


256 


64 


16 
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Table IV (continued) 
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% NORM1 ERROR in IMAGE DOMAIN 


# of radial samples 


1280 


256 


256 


256 


Near. Neigh. 


24.1 


24.0 


24.2 


27.4 


Linear 


19.1 


19.1 


18.4 


21.2 


Cubic B-Spline 


16.1 


16.1 


15.7 


18.9 


Inverse bilinear 


12.5 


11.8 


12.5 


22.3 


FBP 


38.5 


38.5 


38.5 


38.5 



Table V 



25 





% NORM2 ERROR in IMAGE DOMAIN 


# of views 


256 


256 


64 


16 


# of radial samples 


1280 


256 


256 


256 


Near. Neigh. 


93.3 


93.1 


87.5 


108.4 


Linear 


52.8 


52.6 


50.1 


77.6 


Cubic B-Spline 


50.8 


50.7 


49.4 


75.7 


Inverse bilinear 


40.0 


40.3 


45.8 


89.9 


FBP 


84.5 


84.5 


84.5 


84.5 



Table VI 





% NORMoo ERROR in IMAGE DOMAIN 


# of views 


256 


256 


64 


16 | 


# of radial samples 


1280 


256 


256 


256 


Near. Neigh. 


82.0 


81.7 


93.8 


133.5 


Linear 


67.8 


67.5 


78.4 


121.5 


Cubic B-Spline 


66.5 


66.2 


77.7 


117.0 


Inverse bilinear 


61.8 


61.8 


68.0 


89.3 


FBP 


97.8 


97.8 


97.8 


97.8 



[0074] It may be observed that increasing the radial sampling does not reduce the error in the image domain. How- 
ever, decreasing the radial sampling does introduce errors because the rectangular grid spacing are no longer on the 
so same order as radial sampling. However, the value of the "z step" 36 is fixed because of bandwidth constraints. Thus, 
"N" samples are available for each projection such that N multiplied by "z_step" = 2A. 

[0075] In each case the inverse bilinear interpolation has less error both in the 2D FT domain and in the image 
domain compared to the other listed interpolation methods. The number of views in the 2D FT domain and the image 
domain may also be reduced from 256 without significantly increasing the image error as, for example, can be observed 
55 by the data in Table VI. 
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Claims 

1. A method of reconstructing x-ray tomographic images comprising data from a circular orbit and data from a linear 
orbit, based on a shadow zone formed by projecting data of the circular orbit, such method comprising the steps of: 

5 

representing the image obtained from the circular orbit in a circular horn-shaped volume; 

representing the image obtained from the linear orbit in a corresponding linear horn-shaped volume; 

w converting 2D images in the vertical planes of each said circular horn-shaped volume and each said linear 

horn-shaped volume into a plurality of 2D Fourier planes; 

removing projections of a 3D shadow cone in each of said 2D Fourier planes from the circular orbit data and 
replacing said projections with the corresponding projections obtained from the linear orbit data, so as to gen- 
15 erate a plurality of repaired 2D Fourier planes; 

conducting an inverse 2D Fourier Transform on said plurality of repaired 2D Fourier planes to obtain a respec- 
tive plurality of 2D image planes, representing a resulting horn-shaped volume; and 

20 converting said resulting horn-shaped volume back to the x-ray tomographic image. 

2. The method as recited in claim 1 , wherein said step of representing the image obtained from the circular orbit fur- 
ther comprises the step of selecting a shadow cone having a subtended angle of 20 o , wherein said angle 0 O , is 
defined as the angle which identifies the limit of said shadow zone in radon space. 

25 

3. The method as recited in claim 1 , wherein said step of representing the image obtained from the linear orbit further 
comprises the step of generating a plurality of vertical image slices from the linear orbit data to form said linear 
horn-shaped volume. 

30 4. The method as recited in claim 3, wherein said step of generating a plurality of vertical image slices further com- 
prises the step of reconstructing a plurality of parallel beams generated from a projection of the linear orbit data. 

5. The method as recited in claim 4, further comprising the step of using a filtered back-projection method to recon- 
struct said plurality of vertical image slices. 

35 

6. The method as recited in claim 4, further comprising the step of using a Direct Fourier Method to reconstruct said 
plurality of vertical image slices. 

7. The method as recited in claim 6, wherein said step of using a Direct Fourier Method, further comprises the step of 
40 using inverse bilinear interpolation to convert data from the polar grid to the rectangular grid. 

8. The method as recited in claim 7, wherein said step of performing inverse bilinear interpolation further comprises 
the step of performing inverse bilinear interpolation on each line orbit Fourier datum value "x u with respect to the 
four nearest grid locations, gy, gj+1j, gy+1 and g t +1 j + 1 . 

45 

9. The method as recited in claim 8, wherein said step of performing inverse bilinear interpolation on each said line 
orbit datum value, "x," further comprises the step of determining a normalized sub-area, a 00 , a 01 , a 10 , 1t between 
a location of each respective said line orbit datum value M x, u and each of the said four nearest grid locations, gy, 
gi+lj. gy+1 and g,+1j +1 . 

50 

10. The method as recited in claim 9, further comprising the step of collecting for each grid point, g if j, adjacent to a 
respective said line orbit datum value, n x u and assigning a weight, Wy. 

11. The method as recited in claim 10, wherein said step of determining a normalized sub-area further comprises the 
55 step of updating each said grid location according to the following equalities: 

9ij = 9ij + a 11* 9mj= + a oi* 
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12. The method as recited in claim 11, further comprising the step of normalizing each grid data, gy by dividing a 
respective gy by Wy. 

5 

13. An apparatus (10) for reconstructing x-ray tomographic images comprising data from a circular orbit (28) and data 
from a linear orbit, based on a shadow zone formed by projecting data of the circular orbit, such apparatus com- 
prising: 

w means for representing the image obtained from the circular orbit in a circular horn-shaped volume (28); 

means for representing the image obtained from the linear orbit in a corresponding linear horn-shaped volume; 

means for converting 2D images in the vertical planes of each said circular horn-shaped volume and each said 
15 linear horn-shaped volume into a plurality of 2D Fourier planes; 

means for removing projections of a 3D shadow cone in each of said 2D Fourier planes from the circular orbit 
data and replacing said projections with the corresponding projections obtained from the linear orbit data, so 
as to generate a plurality of repaired 2D Fourier planes; 

20 

means for conducting an inverse 2D Fourier Transform on said plurality of repaired 2D Fourier planes to obtain 
a respective plurality of 2D image planes, representing a resulting horn-shaped volume; and 

means for converting said resulting horn-shaped volume back to the x-ray tomographic image (28). 

25 

14. The apparatus (10) as recited in claim 13, wherein said means for representing the image obtained from the circu- 
lar orbit (28) further comprises means for selecting a shadow cone (30) having a subtended angle of 29 0 , wherein 
said angle 0 0> is defined as the angle which identifies the limit of said shadow zone in radon space. 

30 15. The apparatus as recited in claim 13, wherein said means for representing the image obtained from the linear orbit 
further comprises means for generating a plurality of vertical image slices from the linear orbit data to form said lin- 
ear horn-shaped volume. 



16. The apparatus as recited in claim 15, wherein said means for generating a plurality of vertical image slices further 
35 comprises means for reconstructing a plurality of parallel beams generated from a projection of the linear orbit data. 

17. The apparatus as recited in claim 16, further comprising means for using a filtered back-projection method to 
reconstruct said plurality of vertical image slices. 

40 18. The apparatus as recited in claim 16, further comprising means for using a Direct Fourier Method to reconstruct 
said plurality of vertical image slices. 

19. The apparatus as recited in claim 18, wherein said Direct Fourier Method means, further comprises means for 
using inverse bilinear interpolation to convert data from the polar grid to the rectangular grid. 

45 
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FIG. 6 
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FIG. 8 
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